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ABSTRACT 



, ^ ! Symplectic integration algorithms are well-suited for long-term integrations of 

Yf^ , Hamiltonian systems because they preserve the geometric structure of the Hamil- 

^~~^ ', tonian flow. However, this desirable property is generally lost when adaptive 

^ ! timestep control is added to a symplectic integrator. We describe an adaptive- 

^ ! timestep symplectic integrator that can be used if the Hamiltonian is the sum of 

CN '. kinetic and potential energy components and the required timestep depends only 

Fq ! on the potential energy (e.g. test-particle integrations in fixed potentials). In 

O ! particular, we describe an explicit, reversible, symplectic, leapfrog integrator for 

Q>^ ! a test particle in a near-Keplerian potential; this integrator has timestep propor- 

r-| ! tional to distance from the attracting mass and has the remarkable property of 

^i integrating orbits in an inverse-square force field with only "along-track" errors; 

O I i.e. the phase-space shape of a Keplerian orbit is reproduced exactly, but the 

^ ' orbital period is in error by 0(A^~^), where A^ is the number of steps per period, 
cd ■ 

>'■ 

K^ , 1. Introduction 

c3 I During the last decade a great deal of effort has been devoted to the de- 
velopment of symplectic integration algorithms (siAs) for Hamiltonian systems 
( lUhannell fc Scovel 199U| , [Yoshida 1993i |Marsden, Patrick fc Shadwick 19961) . An 



SIA is a symplectic mapping of phase space z = (q, p) and time, M/^ : (z, t) —>■ 
{z',t' = t + h), that approximates the Hamiltonian flow over a small interval 
h. SlAs preserve much of the geometric structure of the Hamiltonian flow; as a 
result, they usually have only oscillatory and not secular errors in the integrals 
of motion and are useful when the main goal is minimizing long-term qualitative 
errors rather than achieving the highest possible short-term precision. 

The most popular SIA is the leapfrog or Verlet method, which can be applied 
to separable Hamiltonians of the form 

i7(q,p,t)=T(p) + t/(q,t) (1) 
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(usually T and U are the kinetic and potential energy). We may define "drift" 
and "kick" operators 

Bh : (q, P, t) ^ ( q + h—, p,t + h\ , K,, : (q, p, t) ^ ( q, p - h—, t J , 

(2) 
and tlie leapfrog operator is then 

L/j = D/j/2K/jD/j/2, or Lh = K/j/sD/^K/j/a; (3) 

either of which defines a second-order SIA ("DKD leapfrog" and "KDK" leapfrog). 
Higher-order SIAs can be derived by concatenating leapfrog operators; for exam- 
ple L^./jL^o/^Lxifc is a fourth-order SIA if a;o = -2^/3/(2 - 2^/^), xi = 1/(2 - 2^/3) 
(|Yoshida 199U|). Leapfrog and its higher-order generalizations have several ap- 



pealing features: only a small number of force evaluations is required per step (1 
for a second-order integrator, 3 for a fourth-order integrator); no auxiliary vari- 
ables are required, thus minimizing memory requirements; and the integrators 
are explicit and time-reversible. 

One limitation of SlAs is that they are usually restricted to a fixed timestep. 
When a standard adaptive-timestep prescription is applied to an SIA, its perfor- 
mance is no better than that of non-symplectic integrators; the reason is that the 
mapping M/j(z)(z,t) is not symplectic even when M/j(z,t) is. This is a serious 
limitation, since in most applications a fixed timestep is inefficient. 

There have been many attempts to construct adaptive-timestep siAs. In the 
context of molecular dynamics, Skeel & Biesiadecki (1994) split the interaction 
potential into a short-range component (rapidly varying, cheap to calculate, zero 
outside a limited range) and a long-range component (slowly varying, expensive 
to calculate); they then evaluate the effects of the short-range potential every 
timestep using a symplectic integrator, adding in N times the long-range po- 
tential at every A^*'^ timestep. This procedure retains symplecticity and can be 
generalized by splitting the potential into any number of parts. Duncan, Levison 
& Lee (1998) describe a variant of the Skeel-Biesiadecki method for long-term 
integrations of solar-system orbits. 

In some situations an adaptive-timestep integrator can be replaced by an 
integrator that uses different — but constant — timesteps for different subsystems. 
In the context of solar system integrations, Saha & Tremaine (1994) describe 
an SIA that uses a different timestep for each planet, which works well so long as 
the planetary orbits are well-separated. 

An alternative to adaptive-timestep SlAs is to abandon symplecticity but de- 
mand that the integrator is time-reversible. Formally, an integrator is reversible if 
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M/jTM/i = T where T is the time-reversal operator. Reversible maps have many 
of the same geometric properties as symplectic maps (Arnold 1984 calls the simi- 
larities "astonishing" ) and hence we expect that reversible integration algorithms 
have similar virtues to SlAs when they are applied to reversible Hamiltonian sys- 
tems. In particular we expect that reversible methods should not exhibit secular 
errors in the integrals of motion. 

Reversible integration algorithms with adaptive timestep are relatively easy to 
generate. Any integration algorithm N/^ can be converted to many reversible algo- 
rithms of the same order ( |Hut et al. 1997|) ; one example is M/^ = N/j/2TN^,^2T- 



Moreover any reversible integration algorithm remains reversible with variable 
timestep if the timestep depends symmetrically on the initial and final phase- 
space coordinates ( [Hut, Makino fc McMillan 1995| ); unfortunately such integra- 



tors are usually implicit and therefore slower than explicit integrators. Various 
explicit reversible integrators with adaptive timestep are described by Huang 
& Leimkuhler (1997), Quinn et al. (1997), Calvo et al. (1998), and Evans & 
Tremaine (1999). 

One problem that requires adaptive timestep is the integration of highly ec- 
centric near-Keplerian orbits (e.g. long-period comets, which often have eccen- 
tricities e > 0.99995). Here the standard approach is to regularize the equations 
of motion; in particular the Kustaanheimo-Stiefel (KS) regularization converts 
the Keplerian Hamiltonian to a harmonic oscillator Hamiltonian, which is easier 
to integrate numerically. KS regularization can be extended to handle few-body 
problems but is restricted to inverse-square interparticle forces. This technique is 
also useful in simulations of star clusters, to deal with the delicate problem of the 
formation and dynamical evolution of tightly-bound binary stars over very long 
times (e.g. [Mikkola fc Aarseth 19931 ). Mikkola (1997, see also [Ranch fc Holman] 



1999[) has combined KS regularization with an efficient SIA designed for nearly 



integrable problems by Wisdom & Holman (1991) to provide a sophisticated 
integrator for eccentric near-Keplerian orbits. 

A closely related problem in numerical celestial mechanics is the long-term 
integration of moderate-eccentricity planet-crossing orbits (e.g. Earth-crossing 
asteroids. Centaurs, Jupiter-family comets), which are nearly Keplerian for mil- 
lions of years, yet occasionally suffer strong perturbations from close planetary 
encounters that may only last a few hours ( Puncan et al. 1998[) . 

The aim of this paper is to discuss a class of explicit adaptive-timestep SlAs 
that can be used to follow Hamiltonian systems of the form (|l]) in the important 
special case where the timestep depends only on the potential energy t/(q, t). We 
provide a brief review of SlAs in §E|, and discuss the use of extended phase space 
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to derive adaptive-timestep SlAs in §^. A class of explicit adaptive-timestep SlAs 
that are particularly suitable for following orbits in nearly Keplerian potentials 
is presented in §H, and numerical tests are described in §^. Section ^ contains a 
brief summary. 

As this work neared completion, we learned of a paper by Mikkola & Tanikawa 
(1999) that contains many similar conclusions. 

2. Review of symplectic integration 

A SIA is a mapping of the form M/^ : (z,t) — » (z',t' = t + h) where M/j 
is symplectic and z' approximates the phase-space trajectory of z over time h, 
which is given by 

i = {z,H}, (4) 

where H{z,t) is the Hamiltonian and the braces stand for the Poisson bracket. 

The SIA can be defined implicitly by a generating function W = W{c{,p',t), 
with the equations of transformation 

dW , dW 

In the simplest SIA, the generating function is chosen to be 

iy = q-p' + /ii7(q,p',t), (6) 

which implies 

9i7(q,p',t) 



q = q + /i- 



ap' 



P = P-^ gq • (7) 

These equations define an implicit first-order SIA. Higher order schemes can be 
derived from more complicated generating functions (e.g. phannell fc Scove] 
TOT. 



We can go further if the Hamiltonian is separable and autonomous, that is, if 

H = Ha + Hb (8) 

where Ha and Hb are time-independent and separately integrable. For a system 
of this type, the equations of motion can be written as 

z = {z,Ha + Hb}. (9) 
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We introduce the differential operators A = {•, Ha}, B = {-, Hb}, and write the 
formal solution of equation (||) as 

z(t) =exp[t(A + B)]z(0). (10) 

By assumption, we know how to calculate exp(tA) and exp(tB). In general, these 
operators are non-commutative so exp[t(A + B)] 7^ exp(tA) exp(tB). The correct 
expression is given by the Baker- Campbell-Hausdorff (BCH) identity ( |Yoshida 



exp(X) exp(Y) = exp(Z), (11) 

where Z = X + Y + i[X, Y] + ^^[X - Y, [X, Y]] H . To construct an explicit 

symplectic integrator of order n we use the BCH identity to find a set of real 
numbers {ci,di) such that 

k 

exp[/i(A + B)] = n expiahA) exp{dihB) + 0(/i"+^); (12) 



the integrator is then 



z(t = h) 



Y[ exp(cj/iA) exp((ij/iB) 



i=l 



z(0), (13) 



where the operators are applied in the order (ci, di, . . . , c^, d^); this is an example 
of the general technique of operator splitting. See Yoshida (1990) for a system- 
atic strategy of finding these numerical coefficients. In the important case of a 
Hamiltonian of the form (|l|), this map takes the simple form (cf. eq. |^) 

-Q- ) , Pi+1 =Vi-hdi{ — \ , (14) 

where z(0) = zi, z(/i) = z^+i. The usual second-order leapfrog (eq. |^) corre- 
sponds to the choice A; = 2, ci = C2 = |, (ii = 1 and ^2 = for DKD leapfrog (or 
ci = 0, C2 = 1, (ii = (^2 = I for KDK leapfrog). Using the BCH identity, one can 
show that 

exp(|/iA) exp(/iB) exp(i/iA) = exp[/i(A+B + C)] where C = {-, i^err} (15) 

and 



12 
Thus DKD leapfrog describes the equations of motion in a surrogate Hamiltonian 



i/err = TT^K^^A, Hb}, Hb + \Ha} + 0(/i^). (16) 



H = H + H,,,; (17) 
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(inore precisely the numerical trajectory lies exponentially close to the exact 
trajectory of the surrogate Hamiltonian) . The good properties of symplectic 
integrators, such as the absence of secular errors in the energy, are a consequence 
of the existence of this Hamiltonian. 



3. Variable timestep and extended phase space 

We want to construct an explicit adaptive-timestep SIA to integrate the Hamil- 
tonian (|l]). Following Mikkola (1997), we extend phase space by introducing a 
fictitious time variable r through the relation 

dt = g{q,p,t)dT (18) 

and take t = Qq as a new coordinate together with the corresponding conjugate 
momentum pq = —H. Thus an extended phase space is defined by 

Q = (go,q), P = (po,p), (19) 

and the equations of motion in the extended phase space are 

dQ , ,dH dV 

dP , ,dH dV 

where 

r(Q, P) = g{q, p, qo)[H{q, p, go) + Po], (21) 

and only trajectories on the hypersurface F = in the extended phase space 
correspond to solutions of the equations of motion in the original phase space. 
The equations of motion (pOD are Hamiltonian in the extended phase space if 
F is chosen as the Hamiltonian. We can now integrate the equations of motion 
with an SIA having constant fictitious timestep At in the extended phase space, 
which is equivalent to a variable timestep At = g{q, p, t) Ar in the reduced phase 
space. 

The Hamiltonian ( plf ) is not generally separable and hence operator-splitting 
techniques cannot be used to derive exphcit SlAs with arbitrary timesteps. Nev- 
ertheless, this approach can yield useful SlAs for specific choices of the timestep 
function g{q, p,t). 
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3.1. A separable Hamiltonian in extended phase space 

We choose the timestep function to be 

.n p^ /(Te(P))-/(-[/(Q)) 

where Te(P) = T(p) + po and ?7(Q) = U{q,qo). The Hamiltonian (|2T|) becomes 

r(Q,P) = /(Te(P))-/(-f/(Q)), (23) 

which is separable. The equations of motion are 

dq^ dV dTip) 

T = ^ = /Wp)+Po), 

(XT Opo 

dr dqi dqi 

|^ = -| = -n-UMf-^. (24) 

To choose the function / we recall that -ff(q, p, go) +Po = 2^e(P) + t^(Q) = 
for the Hamiltonian flow, i.e., Te(P) ^ — f/(Q) during the numerical integration. 
Consequently, /(Tg) — f{—U) ^ and we can Taylor expand the function / 
around Te(P) to obtain 

/(Te(P)) = /(-f/(Q)) + [T,(P) + f/(Q)]/'(-t/(Q)) + 0(T, + uy. (25) 

Therefore, equation ( ^21) yields 

g{Q,P)^f'{-U) (26) 

along the integration path. Thus the timestep can be chosen to be an arbitrary 
function of the potential energy, g = g{—U), and a suitable f{—U) is determined 
by integrating g{—U). 

The choice of the timestep function is crucial to the success of an integrator, 
and the restriction that this function depends only on the potential energy U is 
severe. Nevertheless timestep functions of this form can be useful for a variety 
of dynamical problems. 



3.2. Error analysis 

We now illustrate how to analyze the integration errors that arise when fixed- 
timestep SlAs are used to integrate the equations of motion (p^ ) in extended 
phase space. For simplicity we shall assume that the potential is stationary, 
f/(Q) = f/(q), and restrict our attention to DKD leapfrog. 

The error Hamiltonian for DKD leapfrog applied to the Hamiltonian (PB|) is 



-U^rf [nT,)fp.p, [-n-U)U,U, + f'{-U)U,,] + 0(Ar)^ (27) 

where U^i = dU/dqi and summation over the indices i,j G {1, 2, 3} is assumed. 

Once we have the error Hamiltonian, the numerical error in the energy in 
the original phase space is easy to derive. The integrator accurately follows the 
trajectory of the surrogate Hamiltonian 

f = r + Terr = /(T^) - f{-U) + T,„, (28) 
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thus r is conserved along the numerical trajectory. At the starting point (Qj, P 
r = so r = rerr(Qi, Pi) = Tj throughout the integration. Since T is indepen- 
dent of the coordinate go; the momentum Pq is conserved throughout the integra- 
tion and is therefore equal to minus the initial energy Ei. Thus Tg = T(p) +po = 



AE — f/(q) where AE = E — Ei is the energy error. Equation ( 128|) can now be 
rewritten as 

r, = f{AE - f/(q)) - /(-f/(q)) + rerr(Q, P); (29) 

since the energy error is small we can expand in a Taylor series to obtain 

4. Keplerian two-body problem 

The long-term integration of nearly Keplerian orbits is central to the study of 
solar system dynamics, and the Keplerian two-body problem provides a natural 
laboratory for testing integration algorithms. 

For simplicity we work in two dimensions, and to provide more general for- 
mulae we add an extra potential V^(q, t) to the point-mass potential that defines 
the Keplerian problem. The Hamiltonian for a test particle is thus 



/^ 



H{^,p,t) = lp'-^ + V{q,t), (31) 



d^x 


dv^ 


X 


dV 


<Py 


dVy 


dt^ 


dt 


-- -fJ'^ - 


dx ' 


dt^ 


dt 
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where q = (x, y), p = v = {v^, Vy), r^ = x'^+y'^, and /i is the mass. The equations 
of motion are 

-f^— - 7^- (32) 

r^ oy 

There are two natural choices for the timestep function (^(q, p) when integrat- 
ing bound Keplerian orbits: (i) g oc r^/^ ensures that the timestep is a constant 
fraction of the local free-fall time tg ~ r^/^yU"^/^, so that all phases of highly 
eccentric orbits are followed with the same relative accuracy; (ii) g (x r ensures 
that the coordinate trajectory as a function of the fictitious time is that of a 
harmonic oscillator, so there are no high-frequency harmonics that are difficult 
for numerical integrators to follow. 

We can accommodate both of these choices by assuming that the timestep 
function is a power law in radius, 

g{r) = er^fi'-\ (33) 

where e is a constant that parametrizes the size of the timestep; having introduced 
e we can henceforth set Ar = 1 without loss of generality. Since U ^ —fJ'/r for 
nearly Keplerian orbits, equation (EBI) then suggests that we choose 



liL 



-X 



-<+^ if 7^1, 



f{x) = l 1-7 (34) 

I e/ilogx if 7 = 1. 

The corresponding Hamiltonian is 

^([Te(P)]-^+i-[-f/(Q)]-^+i) if 7 7^1 
r(Q,P) = <| 1-7 (35) 

e/i(log[Te(P)]-log[-f/(Q)]) if 7 = 1. 



For the problem we consider here, f/(Q) = — /i/(x^ + y'^Y^'^ + V{x,y,t) and 

dx 



Te(P) = \{v1 + v'i) +po- The equations of motion (^) in the fictitious time read 



e/i 



dr {\v^^ + \vy^ + pqY 

dy Vy 

-— = e/i 



dr {Ivx"^ + \vy^ + Pq) 

dt 1 

e/i- 



dr ^(1^,2 + 1^^2 + ^^),' 

dvx e/i I jj,x dV' 

dr ~ {jj/r - Vp \ r^ dx 
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dr {ji/r - 1/)T V r3 dy 

dpo e/i dV 



) 



+ 
c)y J 

(36) 



dr {12/r - Vp dt 

These equations can be integrated using leapfrog, since the right-hand sides of 
the first three depend only on momenta in the extended phase space, and the 
latter three on coordinates. For example, if the extra potential V = 0, then DKD 
leapfrog for equations (^) with 7 = 1 can be written 

efiv 

ri/2 = r+ 2 I o ' 

e/i 

tl/2 = ^+ 2 I o ' 

v' + 2po 
e/i ri/2 



v' 



^1/2 



r' = ri/2 + 



t' = ^1/2 + . ,.2'^ o ' (37) 

where r = |r|, and po is an integral of motion equal to minus the initial energy. 
The fictitious time r is equal to the eccentric anomaly u to within a linear trans- 
formation; each step of the integration corresponds to An = e(/i/a)^/^ where 
a = —^fi/E is the semimajor axis. 

More generally, if the attracting mass has a trajectory r^,{t) DKD leapfrog 
with 7 = 1 reads 

e/iv 



ri/2 = r + 

tl/2 = t + 



V — 



f 2 + 2po ' 

f 2 + 2po ' 
e/i[ri/2 - r^(ii/2) 



|ri/2-r^(ti/2)p ' 



Po = Po + 



72 — -1*1^1/2 

efi[ri/2 - r^(ti/2)] • dT^(ti/2)/dt 
|ri/2 -r^(ti/2)P 






K)' + 2p^' 
e/i 

(i;')^2p(," 



^' = ^i/2 + 73T7^^- (3^ 



An appealing feature of equations (|38D is that they contain no square roots 



which are the most time-consuming operation in integrating Keplerian orbits 
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by conventional methods; however, square root evaluations do become necessary 
with this integrator as soon as the non-Keplerian potential V{r, t) is non-zero. 



4.1. Error analysis for the Kepler problem 

To analyze the numerical error of DKD leapfrog with 7 = | we take f{x) 
from equation ( ^4|) and set At = 1 in equation (p7|). The leading term of the 
error Hamiltonian becomes 

rerr(Q,P) '^ ^^ 



24[-t/(q)Te(P)]3 
|2t3/2| VL/|2 - p^p^^-Uf/^U,, - lli-Uy/^ + 3T,i/2](p ■ VUy} . (39) 



For the Kepler problem, U = —fi/r and the error Hamiltonian simplifies to 

(40) 



24T? r r" 



3 ' 2r'^/^ r^/^ 



where Tg = |f ^ + po- 

Similarly, for 7 = 1 the leading term of the error Hamiltonian is 

3 3 

rerr(Q,P) = 24[_^|q^Te(P)]^ [2T,\VU\' + p,p,UU,, - 3(p ■ VUy] ■ (41) 

for the Kepler problem, this simplifies to 



re..(Q,P) = ..,;^,'^° ., . (42) 

12r^(^t;^ +Po) 

This formula leads to a remarkable conclusion, specific to this potential and inte- 
grator. The original phase-space variables (q, p) enter Fgrr in the same combina- 
tion that they enter the Hamiltonian F; in other words the surrogate Hamiltonian 
may be written 

3 
f (Q, P) = F(Q, P) + Fe„(Q, P) = efi log W{Q, P) + ^^^t^^^ p^ , (43) 

where Vr(Q, P) = r(|p^ + po)/fi. Thus the equations of motion (^) for F read 
~ (efi e^HPo \ . . . er e^rpo e^ 

(44) 
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while the equations of motion for F read 

z = {z, r} = — {z, W}, % = — , Po = const. (45) 

Thus the exact trajectory z(r) (Hamihonian F) is the same as the numerical 
trajectory z(r') (Hamihonian F) at a shghtly different fictitious time, where the 
two timescales are related by dr = dr'll — e'^po/{6W'^)] or r = r'(l — e^po) since 
W = 1 on the trajectory. In other words the algorithm (^) follows the Keplerian 
trajectory exactly: the position and velocity are precisely those of the Keplerian 
orbit, and the only error is in the time of arrival at a given location (i.e. the 
only error is "along-track" )Q. Although we have only established this result to 
leading order in the error Hamiltonian, we show in the Appendix that it holds 
at all orders, i.e. for arbitrarily large timesteps. This result also applies in the 
more general case where the attracting mass is in uniform motion rather than 
stationary at the origin (cf. eqs. pSf ). 

We also show in the Appendix that the timing error arising from a step An 
in eccentric anomaly is j^{Au)^/n + 0(Au)^ where n = {fx/a^Y^'^ is the mean 
motion and the error is independent of eccentricity. The fractional error in the 
orbital period is then tt'^/SN'^ where A^ is the number of steps per period (eq. 



4.2. Error analysis for the perturbed Kepler problem 

Because DKD leapfrog with 7 = 1 follows a Keplerian trajectory exactly, 
this method is of particular interest for the perturbed Kepler problem, where the 
extra potential ^(r, t) in equations (|36| ) is small but non-zero. To investigate the 
errors in this case, we take the error Hamiltonian (|4T|), set U = —fJ^/r + V (for 
simplicity we assume that V is stationary), and expand to first order in V: 

2/ipo , 4:poV ,,.12, ,r-W 



3, ,2 



Ferr(Q,P)- '^ 



24(V+po)2 






+ ^^ + 4(y+Po) 



^2 

v^V 3(v-r)V ^v-r „^; 



(46) 



When this is evaluated on the trajectory, we have 



Fe„(Q,P) = -i^eV^ + 



g3r 



12- r- . 24 



- SErV + 4/ir ■ W 



r ViVjVij + rv V 6r(v ■ r)v ■ W 



(47) 



^This result holds only for DKD leapfrog, not KDK leapfrog. 
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The energy error is then (eq. 

AE = ri-rerr(Q,P) _ .^g. 

er 

As r ^ we have v ~ r^^^"^. Thus ii V ^ r'^ as r ^ 0, Fgrr ~ r'^ as well. Then if 
A; > 0, the energy error at close encounters with the attracting mass is 

AE= —, r < u; (49) 

er 

in other words the energy error at close encounters is determined by the initial 
conditions and varies as r~^, independent of the form of the perturbing potential 
at small radii so long as F ^ as r ^ 0. 

The divergence of AE as r — i> 0, even when the perturbing potential V ^ 
as r — i> 0, appears to contradict our proof that the integrator tracks Keplerian 
orbits exactly; the resolution is that the integrator only tracks Keplerian orbits 
on the hypersurface F = in the extended phase space, and numerical errors at 
larger radii perturb the trajectory to the neighboring hypersurface F = Fj. In 
fact it can be shown that in this case the integrator is following a Kepler orbit 
exactly, but for an attracting mass /xexp(Fj/e//). Thus even large energy errors 
at close encounters do not signal a catastrophic failure of the integrator, so long 
as |Fj/e/i| -C 1. 

Moreover there is a simple way to correct these errors. Normally the initial 
value of po is set equal to —E, so that F = 0; instead, we modify the initial value 
of po so that F = F + Ferr is zero. This requires 

Po = -E + ^ [exp(-F,/e/i) - 1] . (50) 

5. Numerical tests 
5.1. Keplerian two-body problem 

We have tested these integration methods by following Keplerian orbits with 
eccentricities e = 0.9, 0.99, 0.999, 0.9999, 0.99999, and 0.999999. Each orbit is 
started at pericenter and followed for 2 x 10^ orbital periods (although a shorter 
integration would have been sufficient, since the energy errors are oscillatory 
rather than growing). We characterize the performance of the integrator by the 
maximum energy error |A£'max/-£'o| = max|(i5 — Eo)/Eo\, as a function of the 
number of steps per orbital period. 

Figure [1| shows the energy error that arises from integrating equations (p^) 
using DKD leapfrog with 7 = | (i-e. timestep oc r'^/^). For comparison we 
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have also shown as open circles the energy error for leapfrog with fixed timestep 
(7 = 0) at eccentricity e = 0.9 and 0.99. Clearly 7 = | provides far more accurate 
integrations than 7 = 0. 

We can compare these energy errors to the analysis of §0!. When evaluated 
on the trajectory [Tg = —U) the error Hamiltonian (|^) is 



3»,3/2 



e r 



24^1/2 



3(r-v) 
2r2 



+ 2E 



e^na? 
24 



ecosM 



,3/2 



Se^ sin^ u 



2fl 



ecosn 



1 1/2 



(51) 

where as usual ra, a, and u are the mean motion, semimajor axis, and eccentric 
anomaly. The energy error is then given by equation (BO), 



AE 



e^n^a^ 



24 



1 



e cos Ui 



3/2 



Se^ sin^ Ui 



ecosu / 2(1 

3e^ sin^ u 



ecos'Uj)^/2(l — ecos'u)^/^ 



2 1 



e cos u 



For high-eccentricity orbits started at pericenter {ui 
\AE\ occurs at m ~ cos~^ e, and is given by 



AE 



E 



16(1 



+ 0[e\l-e) 



(52) 
0), the maximum error 

(53) 



For high-eccentricity orbits started at apocenter (uj = vr), the maximum energy 
error occurs at pericenter. 



AE 



E 



3-2V2(i_e)3/2 



+ 



1-e 



(54) 



These formulae show that pericenter starts lead to smaller energy errors than 
apocenter starts, although in the latter case the errors can be reduced by the use 
of corrected initial conditions (cf. § |4.2|) . 

The number of steps per orbit A^ ~ /q dt/g{r), where g{r) is the timestep 



function 



and P is the orbital period. For 7 



e JO 



df 



K 



^l + ecos/)i/2 e(l + e)V2 Vl + e 



2e 



(55) 



where / is the true anomaly and K is an elliptic integral. Plotting equations 
( p3D and (|55D as a parametric function of e we obtain the dashed lines in Figure 
|l], which agree well with the energy errors from the numerical orbit integrations. 
As we discussed in the previous section, integrating the Keplerian equations of 
motion using DKD leapfrog with 7 = 1 (eq. |37| ) yields even better behavior than 
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Fig. 1. — Maximum fractional energy error over 2 x 10^ orbital periods, as a function of 
number of steps per orbit for the Keplerian two-body problem. The curves correspond 
to eccentricities e = 0.9, 0.99, 0.999, 0.9999, 0.99999, 0.999999. The integrator is standard 
DKD (drift-kick-drift) leapfrog with timestep oc r^/^, following equations ( p3| ) and (|36| ) with 
7 = |. The orbits are started at pericenter. The dashed lines show analytic estimates of 
the energy error from equations (pSf) and (pH). The analogous errors for DKD leapfrog with 
fixed timestep and e = 0.9 and 0.99 are shown as open circles; for larger eccentricities the 
fixed-timestep errors are off-scale. 
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7 = |, as in this case the energy error is zero — in fact there is zero error in all of 
the phase-space functions that are constants of motion in a point-mass potential 
(energy, angular momentum and Runge-Lenz vector). To test the practical value 
of this algorithm we must therefore turn to more general Hamiltonians, which 
we now do. 



5.2. The Stark problem 

The Stark problem is to follow the motion of a test particle subject to an 
inverse- square force plus a constant force; the Stark Hamiltonian is 

ff = ip2 _ /^ _ S . r, (56) 

where the Stark vector S is a constant. The Stark Hamiltonian has three con- 
stants of motion and thus is integrable: these are the energy E, the angular 
momentum component along S, and a third analytic integral that we see no 
point in writing out ( Pars 1965| , [Landau fc Lifshitz 1976|) . We restrict ourselves 



to the planar case, which is particularly challenging because all orbits oscillate 
between retrograde and prograde and hence pass arbitrarily close to the attract- 
ing mass. We shall examine only a single integrator, DKD leapfrog, applied to 
the equations of motion for 7 = 1 (eqs. ^), since in this case the trajectory is 
followed exactly when S = 0. 

Ranch & Holman (1999) have recently tested several integrators on the Stark 
problem, and we shall usually use their initial conditions: the initial eccentricity 
e = 0.9, the Stark vector is oriented 45° to the initial line of apsides, and the 
orbit is started at apocenter. The strength of the Stark perturbation is written 
S = rjE^/^ where E is the energy and r^ -C 1 for nearly Keplerian motion. 

The error Hamiltonian and the expected energy error are given by equations 
( ^71 ) and ( ^8] ) with V = — S • r. Figure ^ verifies the functional form AE oc r~^ 
predicted by (P5| ) and demonstrates the improvement during close encounters 
that results from using the corrected initial condition (pO|). 

Figure ^ shows the energy error as a function of steps per orbit A^ and strength 
of the Stark parameter, for integrations lasting n = 10^ orbital periods using the 
Ranch- Holman initial conditions corrected by equation (|50D . We have plotted 
both the maximum energy error, which is dominated by very close encounters, 
and the average of the absolute value of the energy error, which provides a better 
estimate of the typical error. The average error exhibits the A^~^ dependence 
expected for a second-order integrator; the maximum error is much larger and 
more irregular, reflecting its dependence on rare and rather unphysical close 
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Fig. 2. — Fractional energy error as a function of distance from the attracting mass, for a 
numerical integration of the Stark problem using DKD leapfrog with 7 = 1, which integrates 
Keplerian orbits with zero energy error. The integration lasts for 1000 Keplerian periods of 
the initial orbit and the error is plotted every 100 timesteps. The integration parameters 
are /x = 1, e = 0.1, the initial eccentricity is e = 0.9, the Stark vector is 45° from the initial 
line of apsides, and its magnitude is S" = rjE'^/ij, where 77 = 4 x 10~^. For r <^ 1 the points 
lie on a straight line, consistent with the prediction of equation (|49| ) marked by solid line 
segments. The open circles show the much smaller energy errors when the initial conditions 
are corrected using equation (pH). 
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encounters (the typical maximum eccentricity in an integration of this length is 
given by 1 — Cmax ~ "^"^ = 10~^). We have also plotted one error curve for 
uncorrected initial conditions; we see that the correction reduces the errors by 
about one order of magnitude in this case. 

Ranch & Holman (1999) tested the Wisdom & Holman (1991) integrator on 
the Stark problem with rj = 4 x 10~^. They found that the Wisdom-Holman 
integrator — which works very well for low-eccentricity orbits — was generally un- 
stable, in that the energy error grew by a random walk until the orbit escaped to 
infinity. The instability arose through numerical chaos caused by the overlap of 
resonances in the error Hamiltonian, and could only be evaded if the (constant) 
timestep was in resonance with the orbital period, or is small enough that the 
pericenter passage is well-resolved — even though the Wisdom-Holman integrator 
follows Keplerian orbits exactly for any timestep. Our integrator is evidently not 
subject to these limitations. 

Ranch & Holman also tested several other methods. In particular, Mikkola's 
regularized version of the Wisdom-Holman mapping was completely stable, and 
gave energy errors comparable to those shown for our integrator in Figure Q. 
However, we expect that our method is faster in practice because it requires 
fewer calculations per integration step. 

6. Summary 

We have constructed adaptive-timestep, reversible, explicit SlAs for separable 
Hamiltonians of the form (|I|), using extended phase space ( [Mikkola 1997] ); the 
principal restriction is that the timestep must be a function of the potential 
energy alone. 

Integrators of this kind would require modifications for problems with many 
degrees of freedom since the total potential energy of the system is insensitive 
to local conditions that may demand a short timestep (e.g. a close encounter 
between two bodies). However, for test-particle integrations in fixed, smooth 
potentials or few-body systems with similar masses, these integrators can pro- 
vide both adaptive timestep control and the excellent long-term error control 
associated with symplectic and reversible integration algorithms. 

For close encounters or eccentric orbits in few-body gravitating systems, these 
adaptive-timestep SlAs provide an attractive alternative to fixed-timestep SlAs 
in regularized coordinates; moreover, unlike regularized integrators, adaptive- 
timestep SlAs can also be used to follow orbits in non-Keplerian potentials (e.g. 
galaxies) . 
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Fig. 3. — Fractional energy error for the Stark problem, as a function of number of force 
evaluations per orbit. The initial eccentricity is e = 0.9, the Stark vector is 45° from the 
initial line of apsides, and the orbit starts at apocenter and is followed for 10"^ periods. The 
magnitude of the Stark vector is S" = r]E^/fi where r] = 0.001 (solid lines), 0.005 (dashed 
lines), and 0.02 (dash-dot lines). The lower curves (solid circles) represent the average of the 
absolute value of the energy error, and the upper curves the maximum error. The integrator 
is DKD leapfrog with 7 = 1 and initial conditions corrected using equation (|50|) ; the single 
solid curve with open circles represents the average error for rj = 0.001 if no corrector is 
applied. 
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Although we have discussed only second-order leapfrog integrators, higher- 
order adaptive-timestep SlAs can be derived by concatenating leapfrog steps of 
different lengths ( |Yoshida 199CI| ). 



A particularly interesting example of these integrators is offered by equations 
7|), which follow a Keplerian trajectory exactly, with only "along-track" errors. 



This research was supported in part by NASA Grant NAGS- 73 10. We thank 
Seppo Mikkola for sending us a copy of his paper ([Mikkola fc Tanikawa 1999|) , 
which also shows that the integrator (|3^) is exact for Keplerian orbits. We also 
thank Kevin Ranch for thoughtful comments. 
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A. Appendix: An exact integrator for the Kepler problem 

We prove that the 7 = 1 leapfrog integrator ( pS]) follows a Keplerian orbit 
exactly, except for errors in the time. For simplicity we shall assume that the 
attracting mass is at rest at the origin, as in equations (p?]); the extension to an 
attracting mass in uniform motion is straightforward. 

Let (r, v) and (r', v') be the position and velocity at two points on a bound 
Keplerian orbit with eceentric anomalies u and u' respectively. Then (r', v') 
satisfies the relation 

r' = f{u,u')r + g{u,u')v 

v' = ft{u,u')r + gt{u,u')y^, (Al) 



where 



f{u,u') 



cos(u' — u) — e cos u 



1 — e cos u 



g{u,u') = —[sm{u' — u) — esinu' + esmu] 
n 

nsiniu' — u) 



ftiu,u') = -— , 

(1 — ecos-u'j^l — ecosMj 

,, cos(u' — u) — ecosu' ,. , 

9tiu,u') = ^ '- (A2) 



are Gauss's / and g functions (e.g. Danby 1988 ); here e is the eccentricity, 
n = (yu/a^)^/^ is the mean motion and a is the semimajor axis, which is related 
to the radius through r = a{l — ecosu). 
These equations can be rewritten as 



ri/2 


= r + s(m',m)v. 


v' 


= v + z{u',u)ri/2, 


r' 


= ri/2 - s{u,u')v' 



(A3) 



where 



s{u',u) 



[1 — cos{u' — u)]{l — ecosu) 



nsmiu' — u] 



z(u',u) = )-- -. (A4) 

(1 — ecosu'j(l — ecosMj 
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Comparison to the 7 = 1 leapfrog integrator (^) shows that the two maps 
fr, v) -^ (r\ V) are the same if 



e/i [1 — cos(u' — m)](1 — ecosw) 



p2 + 2po n sin(M' — u) 



e/i nsin[u — u 



q\l2 (1 — e'cosM')(l — ecosu) ' 

e/i [1 — cos('u' — «)](! — ecosw') 



(p')2 -|- 2po nsni{u' — u) 



(A5) 



We use the relations p^ + 2po = 2/i/[a(l — ecosw)], (p')^ + 2po = 2/i/[a(l — ecosw')] 
and square the first of equations ( [A3D to eliminate g^/g. After some algebra we 
find that all of the relations ( [A 51 ) are satisfied if 

„1 — cosAm , , ^, 

e = 2 — — . (A6) 

na sm Am 

where Am = u' — u. Thus, we have proved that our mapping follows the Keplerian 
two-body problem exactly in the original phase space. Although our proof is for 
bound orbits, it is straightforward to show that unbound Keplerian orbits (a < 0) 
are also integrated exactly, with 

cosh Am — 1 
n„a sinh Am ' 

here n„ = {—fi/a^Y^'^ and r = — a(ecoshM — 1). 

We must still establish the relation between the timestep given by the second 
and fifth of equations p8| ) and the actual time At^ required to travel from r to 
r'. The timestep is 

t' = t+At = t+efi I — — —- + -—: — —- ] = t+iea(2+e cos u+e cos u'). (A8) 
'^\p^ + 2po {p'y + 2poJ 2 V 



The relation ( [Aq ) then implies that 

[1 — cos(m' — m)] (2 — e cos m — e cos u') 



nAt 



sin(M' — m) 

1 — cos(m' — m) , 

= 2 ; esmM +esmM. (A9) 

sm(M — m) 

On the other hand Kepler's equation states that the actual timestep is given by 

hAIk = u' — u — e sin u' + e sin u. (AlO) 
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Thus the time error is given by 

1 — POS 7x7/ 

n(At - Atx) = 2^-- Au= M^uf + ^{Auf + 0(Atz)^ (All) 

sm Am 

independent of eccentricity. 

If we take A^ = 27i/Au steps per orbit, the timing error per orbit is 

^^ ""' +0{N-'), (A12) 



P 3N^ 
independent of eccentricity. 
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